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A generic method to investigate many-body continuous-variable systems is pedagogically pre- 
sented. It is based on the notion of matrix product states (so-called MPS) and the algorithms 
thereof. The method is quite versatile and can be applied to a wide variety of situations. As a first 
test, we show how it provides reliable results in the computation of fundamental properties of a 
chain of quantum harmonic oscillators achieving off-critical and critical relative errors of the order 
of 10~ 8 and 10 -4 respectively. Next, we use it to study the ground state properties of the quantum 
rotor model in one spatial dimension, a model that can be mapped to the Mott insulator limit of 
the I-dimensional Bose-Hubbard model. At the quantum critical point, the central charge associ- 
ated to the underlying conformal field theory can be computed with good accuracy by measuring 
the finite-size corrections of the ground state energy. Examples of MPS-computations both in the 
finite-size regime and in the thermodynamic limit are given. The precision of our results are found 
to be comparable to those previously encountered in the MPS studies of, for instance, quantum spin 
chains. Finally, we present a spin-off application: an iterative technique to efficiently get numerical 
solutions of partial differential equations of many variables. We illustrate this technique by solving 
Poisson-like equations with precisions of the order of 10 -7 . 

I. INTRODUCTION AND MAIN RESULTS 

Generally, the behaviour of quantum many-body systems cannot be grasped with a purely analytical approach, 
i.e. merely from the properties of its constituents. As a result, studying them often calls for approximation schemes 
that exhibit two features: (i) to provide an effective description of the system based on the identification of relevant 
degrees of freedom, (ii) to prescribe a (numerical) method to predict mean values of observables from this description. 
One such scheme is that based on the Density Matrix Renormalisation Group (DMRG) [ij, l2| , whose applications are 
countless (see for example Q and references therein for a review). 

Recently, a lot of attention has been devoted to regarding DMRG from a quantum information perspective. The 
aim of this effort is twofold. First, to understand when and why DMRG succeeds in describing a many-body system. 
Second, to improve and generalise the method. This programme has met some success. We now understand better 
the relation between the difficulty to describe the state of a many-body system and the entanglement content of this 
state [1, H, @. It is also clearer why DMRG performs poorly in certain circumstances. The cases of spin chains with 
periodic boundary conditions Q or two-dimensional spin lattices Q are quite explicit in this respect. Many other 
interesting proposals have been made along those lines such as methods to simulate the evolution (in real or imagi nary 
time) of a many-body system, both in the finite-size re gim e and in the thermodynamic limit @, [l(J, [0, [H, [HI, 0, [la |. 
or methods to efficiently compute a partition function [l6| . The common notion underlying all these progress is that 
of matrix product states (MPS) [13, EH O US HH , that is an ansatz that allows to efficiently describe the state of 
the system when it is finitely correlated. 

The present work is devoted to the study of many-body system with continuous degrees of freedom (or continuous- 
variable quantum many-body systems), using the MPS description and the related algorithms. The main motivation 
of our work is to provide a tool that will eventually allow to study (1+1) dimensional quantum field theories (similar 
work in this direction has been already considered in [22|, [23|, US])- The structure of this paper is as follows. 



In section UH we describe the procedure of local truncation of the Hilbert spaces of each subsystem for a generic 



model that allows to use the MPS machinery. 



In section UTTl we review several numerical techniques based on MPS in a detailed and pedagogical way, both in 
the finite regime and in the thermodynamic limit. 



In section ITVl we apply our methods to study a chain of coupled harmonic oscillators. This model is a lattice 
version of a (l-fl)-dimensional Klein-Gordon field. We compute the energy, the von Neumann entropies of 
blocks and correlators of the ground state. These computations allow us to proof our method. 

Next, in scctionfVl we perform a detailed study of the 1-dimensional quantum rotor model. This model is known 
to be interesting in studying the Superfluid-Mott insulator phase transitions, or arrays of Josephson junctions. 
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• In section IVTl we point out that the MPS algorithms may have applications well beyond the problems for which 
they have been initially designed. Namely, we show how they could be used to efficiently compute numerical 
solutions for a wide variety of partial differential equations. 

• Finally, in section fVIII we summarize the conclusions of this work. 

We wish to notice that our proposal for treating continuous-variable systems with MPS is qualitatively different 
from the Gaussian MPS presented in (25l |. 



II. TRUNCATION OF LOCAL HILBERT SPACES IN A GENERIC MODEL 



Let us consider an isolated system of N identical particles living on some domain T> and are characterised, each, 
by a continuous degree of freedom x € D. 

We suppose that the behaviour of this system is exactly described by the Schrodinger equation, with a time 
independent Hamiltonian H cxact . The ground state of this Hamiltonian reads 

|*0,cxact) = J dX\ . . . dXAT*o,cxact(xi, . . . ,Xjv)|xi, . . . ,X N ). (I) 

The problem that we want to analyse is to find an approximation, l^o); of l^o, exact) in the sense that Eq = 
(* |iJ cxact |^ ) is close to Bo.exact = (*o,exaet |#exact | *o,exact)- For this, we will make two approximations. First, 
let us expand 9 (xx, . . . , x^) in orthogonal sets of functions associated with each variable: 

oo 

*0,oxact(xi, . . . ,X N ) = ^ c(si,---,SN)4>iy(xi)-"<t>i*P(x N ). (2) 
si ...sjv = 1 

We wish to consider two approximations of Eq.|2|). The first one consists in an appropriate truncation on the 
dimension of the local Hilbert spaces. To be precise, for each particle k, we truncate the associated basis to a finite 
set of d orthonormal functions: 

*o,oxact(xi, . . .,x N ) = &o,d(xi, ■ ■ ■ ,%n) + Rest d . (3) 

In this paper paper we shall always consider this local dimension d to be the same for all particles. Notice that it is 
possible to make a clever choice of the local <i-dimensional subspaces according to the particularities of the specific 
situation. In fact, once the relevant truncated subspaces are chosen, it is possible to make them depend on variational 
parameters which can be optimally tuned in order to minimize the error of the approximations. Furthermore, as long 
as we restrict our computations to the truncated basis, H cxact can be replaced by the truncated Hamiltonian 

d d 

H = E E <^x 3 ® • - - O ^^iJH^^tl^ <S> - • - <S) 0i^>|^? ».-.<» ^^X^ 3 <».-•«) ^i-V*]- (4) 

ix,...,is=l ji,— ,Jjv=1 

Considering Hamiltonians of the form 

^ N N 

-ffcxact = 2 E Ti + E V ( Xi > X j)> ( 5 ) 
i=l i,j=l 

we observe that H can always be given the standard form 

N 

H = Y,<g> h i 

7 i=l 

where h] is a hermitean dxd matrix acting on the local Hilbert space at site i. This fact might not be obvious for some 
interactions (think of the Coulomb interaction for instance). We prove it as follows: suppose we have a potential term 
between two particles whose matrix elements read V(i 1 i 2 ) i (j 1 j 2 ) = J dxidx24>l 1 (xi)4>* i (x2)V(xi,X2)4>i 2 ( x i) ( i ) j2( x 2)- 
Performing the appropriate singular value decomposition, one gets 

d 2 

1 = 1 
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which is of the desired form. 

The second approximation is to assume that the coefficients c(s\, . . . ,Sn) ar e given by a product of finite-size 
matrices. At this point, it is possible to choose among several prescriptions which indeed involve different numerical 
optimization algorithms. We will review three of them. 

III. MATRIX PRODUCT STATES AND OPTIMISATION TECHNIQUES 

Let us review briefly in this Section some of the main numerical optimization techniques based on MPS. Here we 
shall make a distinction among three different cases, namely, those of periodic boundary conditions, open boundary 
conditions, and the thermodynamic limit. 

A. Periodic boundary conditions 

Let us associate d square matrices A(k,Sk) for each particle k = I...N. These matrices are assumed to have 
dimensions x x X- The periodic boundary condition ansatz consists in assuming that the coefficient c(si, . . . , sn) can 
be written, in terms of them, as 

c(si, ...,s N ) = tr(A(l, si) • • • A(N, s N j). (7) 

It is sometimes convenient to use diagrammatic notations for tensors, where open indices are represented by open 
legs, and contraction of indices is represented by gluing legs. The matrix elements A(k, Sfe) Q(c _ 1 . ctfc at site k and the 
coefficient c(s%, . . . , sn) can then be represented as indicated in FigJTJ 



(a) (b) 




FIG. 1: (Color online) (a) Diagrammatic representation of the matrix element A(k, s k ,)a k _ 1 ,a k at site k. (b) Diagrammatic 
representation of the MPS structure of the coefficient c(si, . . . , sjv) (periodic boundary conditions). 

The representation given in Eq. ([7]) is efficient in the sense that it involves a relatively small number of independent 
parameters [7]. Indeed, 0(d N ) independent parameters are necessary in order to describe such a general state as 
Eq.([3]). In contrast, the MPS description requires only 0(Nd\ 2 ) parameters, which for fixed x is linear in N. The 
approximation ([7]) can be made arbitrarily precise upon taking x sufficiently large. Actually taking x — d\ N /^ 
allows to describe any state of N particles with local Hilbert spaces of dimension d 7j within this prescription. But 
remarkably, for many Hamiltonians of physical interest, | XI/ o, exact) can be approximated accurately with quite a modest 
value of % [6[. To be precise, x can be related to the entanglement entropy S of a subsystem in such a way that, for the 
case considered here, x > d s ^ 2 . Hence, MPS provide a faithful representation of the ground state of slightly-correlated 
quantum mechanical systems, as is the case of non-critical 1-dimensional quantum spin chains [5|. 

In addition to providing an economical description of quantum states, MPS also allow for an efficient computation 
of expectation values of observables. Indeed suppose that we want to compute the mean value of a tensor product of 
local observables (Oi (g> • • • (g) On}- One directly checks that the following identity holds 



(Oi <E> • • • <8> N ) =tr(Oi---Ojv), 



(8) 
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where Ok — Yl s ' Sk =i( s 'k\Ok\sk)A(k, s' k ) ® A(k, Sk) is a transfer matrix for the local operator Ok at site k (A denotes 
the conjugate matrix of A). These quantities are diagrammatically represented in Figj2] and it is easy to see that in 
this case the computation of the mean value of a tensor product of local operators requires of a computational time 
which grows as 0(Nd\ 5 )- 



(a) (b) 




FIG. 2: (a) The diagrammatic representation of the transfer matrix Ok at site k. (b) The diagrammatic representation of the 
mean value (Oi ® . . . <g> On) = tr(Oi . . . On), where the names of the indices are not explicitly indicated for simplicity. 

In order to compute an approximation to the ground state of the desired Hamiltonian operator H, the variational 
improvement now proceeds as follows. First, one sets all matrices A to some initial values. Then, one successively 
improves the matrices A related to the particle 1, then those of particle 2, . . ., then those of particle N, and then again. 
Sweeping sufficiently many times through the set {1, . . . , N}, one converges to a stationary value of the approximated 
ground state energy E - The improvement of the matrices related to particle k proceeds as follows. One computes 
the matrices TC^ (the effective Hamiltonian matrix) and Af( k > (the normalisation matrix) such that 

<*> = E E E ^M^^^& t _ lA) , (stai _ lttl) i(M fc ) ai _ 1AI o) 

s' k -s k =l /3fc-i,fts=l ajs_i,ak=l 
d X X 

(*o|*o> = E E E ^( fc ! 4)ft- 1 A<^ fc _ 1&) , (Sfcaii _ iafc) A(A ; , Sfc ) afc _ 1 , ajb . (10) 

s' k -s k = l /3fc-i,fts=l a k _ 1 ,a k = l 

The computation of both Ti^ and Af( k > follows the same rules as the computation of expected values in terms 
of transfer matrices, and is straightforward by using the tensor network diagrammatic representation from FigfJ] 
removing from the diagram the appropriate matrices at site k, given that H has been brought to a normal form ([6]). 
Extremising (if) with the constraint (^ol^o) = 1 with respect to the coefficients A(k, Sk)a k _ 1 ,a k leads to having to 
solve the (generalised) eigenvalue problem 

d X 

E E [ H Kft_ 1 &),( St « t - 1 « t ) - M '(4A_ 1 ft),( Sl a k _ 1 a t )] A (' : . S ^n,^ = 0- V4, A_ b ft, (11) 
Sfc=l a fc _i,a fc = l 

where a Lagrange multiplier A has been introduced. The value of this Lagrange multiplier is actually the value of the 
ground state energy, as computed at the considered step of the algorithm. 



B. Open boundary conditions 

We assume again an ansatz of the form ([7]) but where the matrices now have different sizes, depending on the site 
to which they are related. The matrices A(l,s) have size 1 x x, the matrices A(k,s),k = 2, ...,N — 1 have size 
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X x X, and the matrices A(N, s) have size x x 1- This matrix product structure can again be represented in terms of 
a diagrammatic tensor network, as shown in Fig[3J 




FIG. 3: (color online) Diagrammatic representation of the MPS structure of the coefficient c(si, . . . , sn) in terms of open 
boundary conditions. 

As in the case of periodic boundary conditions, this MPS description is efficient as compared to the exact description 
of the d N coefficients c(s±, . . . , sn). Again, we need of 0(Ndx 2 ) parameters to describe the appropriate product of 
matrices, which can be made as precise as desired by increasing the parameter x- I n this case, it is easy to see by 
considering the consecutive Schmidt decompositions of the system that if % = d^/ 2 ] then any state of N particles of 
local dimensionality d can be represented. 

Working with open boundary conditions is a stronger assumption than working with periodic ones. It is also a less 
physical approximation when studying systems described by a translationally invariant Hamiltonian. However, they 
allow to re-parametrise the MPS in a manner that considerably lowers the computational cost of each step of the 
optimisation procedure. This can be achieved thanks to the local gauge freedom existing in the definition 0. Let 
Tfc denote a x x X hivertible matrix. The transformation A(k, s) — > A(k, s)T k , A(k + 1, s) — > Tj^ A(k + 1, s) (k < N) 
leaves invariant the state described by the MPS. 

We start with the site 1 and perform the singular value decomposition: 

min(ti,x) 

A(l, Sl ) hai = £ U^Ay^l- (12) 

where by definition of singular value decompostion E^ 1 ) is a diagonal min(d, x) x min(<i, x) matrix, and and yW 
satisfy U^U^ 1 ' = I m m(d, x ) an d V^V^ = I m in(d,x)- We then perform a gauge transformation 

A(l,s 1 ) l!ai ^U(l,s 1 ) 1 , x = U^ x , 



A(2,s 2 ) aua2 ^B(2,s 2 ) x , a2 = £4!W,M 2 > S 2W 

f-i—l u—l 

Observe that the dimensions of the matrices related to sites 1 and 2 may change along this transformation. Then we 
perform the further singular value decomposition for B{2, s%) 

min(d 2 ,x) 

5(2, S2 ) A , Q2 = £ U^T^V^, 

and make the gauge transformation 
B(2,s 2 ) x , a2 ^U(2,s 2 )^ = U^ x . , 



min(d 2 ,x) X 

A(3, S 3 ) Q2!Q3 -> 5(3,33)^^3 = ^W V v?2 3 A (3> S 3)a 3 ,a 3 - 

V— 1 Ck.2 — 1 

Iterating this procedure up to the site k — 1, the state is represented through the product of matrices 



(7(1, sx) ■ ■ ■ U{k - 1, s k ^)B{k, s k )A(k + 1, s k+1 ) ■ ■ ■ A(N, s N ). 



(13) 
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Then, one applies the same procedure, starting from site N, up to site fe + 1 this time, keeping at each step the matrix 
at the right of the singular decomposition. Eventually, our state is described by a product of matrices of the form 

17(1, «!)•■■ U(k - 1, s k ^)A{k, s k )V(k + 1, s k+1 ) ■ ■ ■ V(N, s N ). (14) 

The principal virtue of this gauge is that (^ol^o) now assumes a simple form, namely, 

d 

<*o|*o) = ]T tv(A*(k,s k )A(k,s k )), (15) 

so that the corresponding normalization matrix is now equal to the identity. This specific property of the chosen 
MPS representation for open boundary conditions allows a more straightforward numerical analysis of the correspond- 
ing variational optimization algorithm to compute good approximations of the ground state of Hamiltonians. The 
variational improvement proceeds now by sweeping back and forth over the N sets of matrices of the MPS optimizing 
over a particular site at each step. However, when optimizing over the set of matrices at site k, we shall consider 
now the MPS representation for open boundary conditions in terms of matrices U(l, si) for sites I = 1, 2, . . . , k — 1, 
V(l,si) for sites I = k + 1, k + 2, . . . , N, and A(k,s k ). The corresponding generalised eigenvalue problem is then 
transformed into a simple eigenvalue problem. Therefore, computing the ground state of the corresponding effective 
Hamiltonian TiP^ directly gives the variational ground state energy when optimizing over the site fc, together with the 
corresponding well-normalized optimal matrix for that site. This can be very efficiently computed by means of large 
sparse- matrix techniques. Moving from site k to site fe+1 then proceeds by means of evaluating the new representation 
of the MPS in terms of U(l, si) for sites I = 1, 2, . . . , k, V(l, s t ) for sites I = k + 2, k + 2, . . . , N, and A(k + 1, s k +i), 
which only involves performing the appropriate singular value decomposition over the just-obtained optimal matrix 
at site k. As a matter of fact, the computation of TL^ at each step (as well as the computation of expected values of 
local observables) is also very much simplified thanks to the considered structure of the open-boundary MPS, taking 
into account the normalization conditions from the left and from the right of the matrices forming the MPS. 

Finally, it is easy to see that considering the above procedures, together with efficient clever contractions of the dif- 
ferent tensor networks that appear in the calculations, the computational running time of the variational optimization 
algorithm over MPS with open boundary conditions grows like 0(N 2 d\ 3 )- 

Let us now present an alternative derivation of the MPS structure for open boundary conditions If we perform 
the Schmidt decomposition between the local system 1 and the remaining N — 1, we can write the state as 

min(rf,x) 

|*o} = £ A(l) Q1 |rW)|4T" n) ) , (16) 

Ql- 1 

where A(l) Ql are the Schmidt coefficients, \rai) and |r<2'" n ^) are the corresponding left and right Schmidt vectors. 
Expressing the left Schmidt vector in terms of the original local basis for system 1, |Wo) can then be written as 

d min(d,x) 

l*o>=£ E r(l, Sl ) lai A(l) Q1 |^ 1 1 ))|r^-")), (17) 

Si — 1 Ql- 1 

where T(l, Si)i ai is the appropriate coefficients of the change of basis. That is, \tc$) — J2 Sl ^0-i s i)iai )• At this 
point, we expand each Schmidt vector \t^[ in the original local basis for system 2, that is, 



\rf n) ) - 



E i*2 > >i"&; B) >- as) 



82=1 



We now write the unnormalised quantum state |u>£ 3 lS2 ™' ) ) in terms of the at most d 2 eigenvectors of the joint reduced 

density matrix for systems (3, . . . , n), that is, in terms of the right Schmidt vectors \tq 2 ) of the particular bipartition 
between the first two local systems and the rest, together with the corresponding Schmidt coefficients A(2) Q2 : 



min(d ,\) 

i3 ;sf)= E r(2, S2 ) QlQ2 A(2) a2 |r( 3 2 -)). (19) 

a 2 =l 
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Replacing the last two expressions into Ea. (fT7| we get 

d min(cZ,x) min(d 2 ,x) 

l*o)= E E E r(l )Sl ) lai A(l) ai r(2, S2 ) aia2 A(2) a J^)®^))|r( 3 2 --»)). (20) 

Sl,S2 = l ai = l a 2 = l 

Iterating the above procedure, we finally get a representation of the quantum state in terms of some tensors T and 
some vectors A: 

l*o,*> = E E r ( 1 > «i)i Q1 A(i) Q1 r(2, S2 ) Q1Q2 A(2) Q2 . . . A(iv - lW^rXiv, s N ) aN _ a \<pW<p<g ■ ■ ■ C). (21) 

{s} {a} 

where the whole set of indices {s} and {a} are added up to their respective allowed values. 

The above decomposition immediately provides the Schmidt vectors A of all the possible contiguous bipartitions of 
the system. In fact, Eq. (|2ip is indeed a reparametrization of an open boundary condition MPS, in a gauge where 

A(k,Sk) ak _ iak = r(fc, si)a k _ iak A(fc) Qfc ■ (22) 

This reparameterization of the MPS is diagrammatically represented in Fig0J We also see that the maximum 




FIG. 4: (color online) The diagrammatic representation of the decomposition of A(k, Sfc) into a tensor T(k,Sk) and a Schmidt 
coefficient A(fc). 

allowed rank of the different indices otk, k — 1, ... ,N — 1, is site-dependent, since the size of the Hilbert spaces 
considered when performing the consecutive Schmidt decompositions depends on the site. In particular, we have that, 
at most, a k = 1, 2, . . . , d k for k = 0, 1, . . . , [N/2\ , and a k = 1, 2, . . . , S N ^ for k = N,N - 1, . . . , [N/2\ + 1. The 
value of the true rank of the indices in the MPS comes dictated by the minimum between these quantities and the 
parameter \- 

In practice, when studying physically relevant states, many of the Schmidt coefficients for the different contiguous 
bipartitions of the system shall be (almost) equal to zero, depending on the particular state being considered. One 
is then interested in bounding the range of the indices by a number \, which now is clearly understood as a measure 
of the maximum contiguous bipartite entanglement in the system, since it is the maximum allowed Schmidt number 
for a contiguous bipartition in our MPS decomposition. Notice that this truncation in the Schmidt ranks imposes 
a strict restriction on the maximal amount of entanglement entropy that the ansatz can handle. To be precise, if 
S = — tr(plog 2/ o) is the von Neumann entropy corresponding to the reduced density matrix describing one of the 
subsystems of a bipartition of the quantum state into two contiguous pieces, then we see that S < log 2 X- This 
property makes these states to be finitely- correlated. 



C. The thermodynamic limit for translationally invariant systems 

In the case of a translationally invariant system, the techniques described in the foregoin g an alysis can be extended 
so as to consider an infinite number of particles, which allows to get rid of finite-size effects [13| . We now assume that 
all the properties of the state are defined in terms of a single set of matrices A(s) = r(s)A which now do not depend 
on the position. 
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It is indeed possible to perform the evolution of such an MPS in the thermodynamic limit both in real and euclidean 
time as driven by a translationally- invariant Hamitonian, as originally explained in 13]. To this end, let us assume 
now that our Hamiltonian is a translationally invariant infinite sum of interaction terms which only involve nearest 
neighbour interactions, 

H = Y / h {Kk+1) - (23) 
k 

The evolution operator e~ lHt is now decomposed by means of a Suzuki- Trotter decomposition [26] ] in terms of unitary 
gates acting only on two sites: 

jj(k,k+l) = e -ih<- k - k + 1 '> st ^ j- 24 s 

where St <^ I, and which have to be applied t/St times. Arranging the above gates into two mutually non-conmuting 
pieces 

U(AB) = ® k U {2k > 2k+1) 

U(BA) = ® k U {2k ^ 2k \ (25) 

we can compute the evolution in the MPS by slightly breaking the translational invariance as follows: since the set 
of gates in U(AB) can be applied at once (the same holds for those in U(BA)), we consider an MPS invariant under 
shifts of two lattice sites, so that 

r(2k,s 2k ) = T(A,s A ) 

T(2k + 1, s 2fe+ i) = T{B,s B ) (26) 
and 

X(2k) = X(A) 

\{2k + l) = X(B), (27) 

for all the possible values of k. The updating of the MPS along the time evolution is then dictated by sequentially 
updating matrices for the A and B sites according to the application of the gates in Eq. p5)l . As explained in [13j ]. 
the updating rule for the matrices of the MPS once a two-particle unitary gate is applied can be done by means of 
the appropriate singular value decomposition followed by a truncation up to x of the rank of the common index that 
joins the matrices at the two sites. 

Let us be more precise with the above statement. In the situation in which a two-body gate J7^' B - ) from is 
applied over two consecutive sites A and B, the updating procedure starts with the following tensor contraction: 

d d 

Q(s' A , 4) = E E u(A,B \ s 'a, s' b ;s a , s B )X(B)T(A, s A )X(A)T(B, s B )X(B). (28) 

S A = 1 S B =1 

By performing the singular value decomposition of the above tensor with respect to the left and right indices, we get 

Q(s' A , s' B ) = A(s' A )X'(A)B(s' B ) = X(B)T'(A, s' A )X'(A)T'(B, s' B )X(B), (29) 

where T'(A, s' A ) = X~ 1 (B)A(s' A ) and T'(B, s' B ) = B(s' B )X~ 1 (B) respectively correspond to the updated new matrices 
for sites A and B. The last step in the updating procedure comes dictated by a truncation in the eigenvalues A' (^4), 
in such a way that the final size of the updated matrices is again \ x X- The sequential application of gates and the 
truncation scheme are diagramatically represented in FigfS] 

Following this updating rule, it is possible to compute the evolution in time of any MPS in the thermodynamic 
limit as driven by a Hamiltonian such as the one from Eq. (|23p in 0(d 3 x 3 ) time. Also, by performing evolution in 
euclidean time and taking very good care of the normalization of the wave function at each computational step, it 
is possible to find MPS approximations to the ground state of a wide variety of interacting Hamiltonians such as 
the ones considered in this paper. Note the difference between this optimization procedure, based on euclidean time 
evolution in the thermodynamic limit, and the ones described for finite systems in terms of local optimizations of 
the matrices by performing sweeps over the different particles, which involved the solution to (generalised) eigenvalue 
problems. 
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(a) 



A B A B A B 




KB) MB)' 1 X{BY l X{B) 



FIG. 5: (color online) (a) The sequential application of two-particle unitary gates to the MPS structure, (b) The updating rule 
for the MPS in the thermodynamic limit under the action of a two-particle gate follows the singular value decomposition (SVD) 
of the composed tensor, and a local truncation up to \ of the common index between the two particles. After the singular value 
decomposition, the vectors A are properly inserted by hand on the two sides in order to preserve the MPS structure, which also 
changes the values of the tensors at A and B by a factor of A -1 . 

IV. BENCHMARK: THE HARMONIC CHAIN 

As an illustration of the above methods, we consider an open chain of spin— particles interacting via the Hamil- 
tonian 

N N N-l 

^oxact = -^p'+M^^+A^^-^+i) 2 , (30) 
i— 1 i— 1 i— 1 

corresponding to the well-known case of interacting harmonic systems [H, HE H3, HE HH 1 aim which is critical 
in the thermodynamic limit for ji = 0. A remarkable feature of this Hamiltonian is obtained upon considering the 
continuous limit 

lim lim ff mct « / ' dq{d v d u + m 2 )ip(q), (31) 

N^ooA— >oo J 

where we have made the identification i — > q, x — > m 2 , and which corresponds to the usual Klein-Gordon 
Hamiltonian of a free spinless field ip. 

The ground state of H cxact and its associated energy, |\&o, exact) and i^o, exact, can be analytically computed [3l| . 
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Rewriting H oxact as 

N N 

■ffcxact = X 51 Pi + X! X * K ij X h ( 32 ) 
»=1 i,j=l 

we have ^o.exactOcij • • • j af/v) = Ce~ Sjj^i 1 '^)^^^ (7 being a normalization constant, and i?o, exact = |tr(\/ir). 
This result will allow us to estimate the accuracy of our calculations. 

Let us now apply the methods discussed in the previous sections and search for an MPS approximation to the 

ground state. As basis functions, we have chosen to use, for all particles, the (truncated) spectrum of a single 
harmonic oscillator. That is, 

3 -(x) = Cj e - a2x2 / 2 Hj{ax), j = 0...d-l. (33) 

In this expression, Cj is a normalisation constant, TCj denotes the degree- j Hermite polynomial, and a = (2/x) 1//4 
(notice that we could also keep a as an external variational parameter). With this choice H cxac t can be approximated 
by the truncated Hamiltonian 



N N-1 N-1 

H = ~J2 r (k) + ( A + ^)(C (1) + C (A °) + (2A + /*) £ C (fc) - 2A ]T 77 (fe V fe+1) , (34) 



k=l k=2 k=l 

where r, £ and rj are hermitean matrices defined as 



'2 



d 

T m , n = j dx (j) m (x)—^(j> n (x) 



dx 4> m (x) X 2 (j) n (x) 

dx4>^{x)x<t> n {x), (35) 



where m, n = . . . d — 1. Analytical expressions for the entries of the matrices r, £ and rj are straightforwardly derived 
from the elementary properties of Hermite polynomials (see [13]). The Hamiltonian (f3"3"|) is now of the form ©, so 
that its properties can be investigated with the algorithms of the previous sections. To be precise, we show here the 
performance of the finite-size numerical optimization algorithms for a system of = 30 sites and open boundary 
conditions, for several values of the external parameters. Notice that our simulations allow to change both d and x 
independently. 

A. Relative error of the ground state energy 

We have performed computations of an approximated ground state of the system and compared the obtained results 
with the exact ones, for several values of d and \- O ur comparison is made precise by means of calculating the relative 
error of the ground state energy, as measured by the quantity 



5 = 



En — El 



0. exact 



Eo,e 



(36) 



In the off-crtical regime (/x 7^ 0), it is seen in table U that the relative error of the ground state energy decreases 
exponentially with d down to 10 -8 , while it remains very stable for the considered values of X- O n the other hand, 
the exponential behaviour of the relative error as a function of d becomes weaker when close to the criticality, as can 
be observed also in table U where it decreases up to 10~ 4 . Let us remark that the precisions that we find in our 
calculations are in agreement and can compete with those found in previous numerical analysis of the harmonic chain 

B. Off-critical correlator 

Our numerical calculations allow us to compare the performance of the simulation technique by means of considering 
the correlation functions of the system. To be precise, we have numerically computed the two-point correlation function 

C{L) = {xlxl)-{xl){xl), (37) 
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N 


d 


X 


6{ji = 1) 


S{fl = 10" 


6 ) 


30 


4 


5 


8.27 x 10" 


3 


5.38 x 10" 


-3 


30 


4 


10 


8.27 x 10" 


3 


5.23 x 10 


-3 


30 


4 


15 


8.27 x 10- 


3 


5.23 x 10 


-3 


30 


4 


20 


8.27 x 10" 


3 


5.23 x 10 


-3 


30 


4 


25 


8.27 x 10" 


3 


5.23 x 10" 


-3 


30 


6 


5 


5.95 x 10" 


4 


3.07 x 10" 


-3 


30 


6 


10 


5.95 x 10" 


4 


1.71 x 10 


-3 


30 


6 


15 


5.95 x 10" 


4 


1.68 x 10 


-3 


30 


6 


20 


5.95 x 10" 


4 


1.67 x 10 


-3 


30 


6 


25 


5.95 x 10" 


4 


1.67 x 10" 


-3 


30 


8 


5 


4.03 x 10" 


5 


3.01 x 10" 


-3 


30 


8 


10 


4.01 x 10" 


5 


1.04 x 10 


-3 


30 


8 


15 


4.01 x 10" 


5 


8.65 x 10 


-4 


30 


8 


20 


4.01 x 10" 


5 


8.48 x 10 


-4 


30 


8 


25 


4.01 x 10" 


5 


8.46 x 10" 


-4 


30 


10 


5 


2.82 x 10" 


6 


2.94 x 10" 


-3 


30 


10 


10 


2.63 x 10" 


6 


9.36 x 10 


-4 


30 


10 


15 


2.63 x 10" 


6 


6.46 x 10 


-4 


30 


10 


20 


2.63 x 10" 


6 


5.34 x 10 


-4 


30 


10 


25 


2.63 x 10" 


6 


5.24 x 10" 


-4 


30 


12 


5 


3.63 x 10" 


7 


3.06 x 10" 


-3 


30 


12 


10 


1.73 x 10" 


7 


9.47 x 10 


-4 


30 


12 


15 


1.72 x 10" 


7 


5.72 x 10 


-4 


30 


12 


20 


1.72 x 10- 


7 


4.4 x 10" 


4 


30 


12 


25 


1.72 x 10" 


7 


3.70 x 10" 


-4 


30 


14 


5 


2.05 x 10" 


7 


2.98 x 10" 


-3 


30 


14 


10 


1.42 x 10" 


8 


9.07 x 10 


-4 


30 


14 


15 


1.41 x 10" 


8 


5.60 x 10 


-4 


30 


14 


20 


1.41 x 10" 


8 


3.75 x 10 


-4 


30 


14 


25 


1.41 x 10" 


8 


3.21 x 10" 


-4 



TABLE I: Performances of the MPS algorithm for computing the ground state energy of the harmonic chain, both in the 
off-critical regime (fi = 1) and close to criticality (// = 10 -6 ). iV denotes the number of variables, d denotes the number of local 
levels, x is t ne dimension of the matrices of the MPS ansatz, and 5 the relative error. 

for N — 30, x — 30 and d = 5, and compared it to its exact value. It is a fact that, when close to criticality, the 
correlators are not well reproduced by means of our MPS techniques (which inherently make use of a finite correlation 
length). However, away from criticality (/i = 1), our algorithms provide good approximations of the correlation 
functions, as can be seen in Fig. [6l where it is shown that the numerical results obtained for C(L) agree very well 
with the exact value. 



C. Absolute error of the entanglement entropy 

Recently, a lot of attention has been devoted to computing the entanglement present in Gaussian systems (see for 
instance [27} and references therein). This motivated us to look at the entanglement entropy S — — tr(plog 2 p) of the 
ground state reduced density matrix p of half the open chain, as a further test of our method. We have computed S 
for x = 10) d = 14 and A = 0.5, as a function of the parameter p 2 , and compared it to the exact case. As can be 
seen in FigJTl the absolute error decreases superexponentially when departing from the critical region, easily achieving 
values of the order of 10 -9 when being sufficiently away from criticality, and giving small relative deviations from the 
exact entropy in FigUl 
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FIG. 6: Off-critical correlator C(L) = (x\ x 2 L ) — {x\)(x 2 L ) for N = 30 and /i — 1. The numerical computations have been done 
for N — 30, x = 30, and d = 5. Very good agreement with the exact value is observed to appear. 

V. THE 1-DIMENSIONAL QUANTUM ROTOR 

As a further example of the application of the techniques presented in the previous sections, we consider here the 
properties of the ground state of the so-called quantum rotor model in 1 spatial dimension. This continuous-variable 
model is defined by means of the Hamiltonian 



N 

£ 

k=l 



-2 J cos (0k — 



'k+U - 77 



UjP_ 

2 de 2 k 



(38) 



where the angles 9k € [0, 27r) Vfc, and U, J are external parameters. The above Hamiltonian accurately describes 
the behaviour within the Mott-insulator regime of the 1-dimensional Bose-Hubbard model, where small fluctuations 
around the mean number of atoms per site can be considered [33l [34[ , hence providing an approximate description 
in a specific regime of the behavior of cold atoms in 1-dimensional optical lattices. This model has also met success 
in reproducing the different properties regarding the tunneling of Cooper pairs between different superconducting 
islands in arrays of Josephson junctions [351 ]. From a technical point of view, note that the model is invariant under 
global rotations of all the angles, that is, under the global U(l) symmetry group. A suitable tuning of the parameters 
defining the model describes then in the thermodynamic limit the evolution towards a quantum phase transition, 
whose {/(l)-invariant critical point can be described in terms of the conformal field theory of a free scalar field with 
central charge c = 1, as proven in [36| . 

Let us perform now a suitable truncation on the local dimension of the different Hilbert spaces as follows. We 
choose the individual particle basis of plane waves 



(39) 



with Sk = —(d — l)/2, . . . , 0, . . . , (d — l)/2, for all the possible k. Restricted to this truncated basis, the Hamiltonian 
from Eq. ({3"8")) can be expressed as 



N 



k=l 



H = J2(- 2J [v {k) V {k+1) + W {k) W {k+1) ) + -T« 



(40) 
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FIG. 7: Absolute error in the entanglement entropy for the harmonic chain for N = 30, d — 14, x — 10 and A = 0.5, as a 
function of fi 2 . 



The matrix elements of operators V, W and T adopt the simple form 



* / m . n. — 



w„ 



lis 



m— l,n 



m— l,n 



Tn 



(41) 



Note that for d = 2 we have that V = ha x , W = \a y and T = (<r z ) 2 . The truncated Hamiltonian given in Eq. (f40|) 
for d = 2 then represents a very specific case of the spin-1/2 XY quantum spin chain, which was already known to 
be an approximation of the 1-dimensional Bose-Hubbard model in the Mott-insulator regime [j| . 

We have performed a variety of calculations over the above truncated Hamiltonian for several values of the local 
dimension d of the Hilbert space, by means of the MPS-like algorithms described at the beginning of this work. Our 
results extend those computed in [33| . and can be understood in terms of the underlying quantum phase transition 
that the model undergoes in the thermodynamic limit. 



A. Ground state energy and entanglement entropy of half a system 



We have performed computations of the ground state energy per lattice site of the system eo, both in the finite-size 
regime and in the thermodynamic limit, as a function of the ratio U / J. Initially, we have performed our analysis for 
a local dimension d = 7 and \ = 20. The results can be found in FigGj] and Figfll]! which agree with the previous 
results from [33| , and where it is possible to see the convergence with N of the ground state energy found with the 
algorithms for MPS with open boundary conditions and for MPS in the thermodynamic limit. 

In order to have a more accurate description of the properties of the ground state close to the critical point, we 
have performed as well a computation of the entanglement entropy of half a system, as a function of the ratio U / J, 
and again both in the finite-size case and the thermodynamic limit. Our results are represented in Fig llll which 
were obtained by using the same values in our MPS techniques than those used in the previous computation of the 
ground-state energy. According to the maximum observed value of the entanglement entropy, the quantum phase 
transition appears to be in our infinite-size simulations around U/J ~ 3, also in accordance with previous estimates 
from 33]. For finite N, our MPS numerical algorithms constantly exhibit spontaneous symmetry breaking of the 
Z2 symmetry of the ground state when approaching the limit U/J — > 0, which makes the entropy suddenly decrease 
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FIG. 8: Exact entanglement entropy for the harmonic chain for N = 30 and A = 0.5, as a function of . 



towards (separable state). The MPS algorithm in the thermodynamic limit has shown to be much more stable in 
this respect, as it always maintains the GHZ-like superposition in the ground state when approaching the same limit, 
implying the entropy to evolve towards 1. A divergence in the correlations is also seen to appear when increasing N 
at the critical point. However, we can only hint such a divergence because of the finite value of x that we use in our 
simulations. 



It is possible to make a check of consistency of the universality class of the model by studying the finite-size 
corrections of the ground state energy when close to the critical point. As proven in [371 ] . for a system with finite size 
iV and open boundary conditions, the ground state energy approaches its critical value in the thermodynamic limit 
in the way dictated by the relation 



where £o(o°) is the energy per particle in the continuum limit, a is a constant, and c is the central charge describing 
the underlying universality class of the corresponding conformal held theory in the thermodynamic limit. Fitting 
the above expression to the ground-state energies computed at the value of maximum entropy for finite N can then 
provide us with an estimation of the value for the corresponding central charge c. Our results for the 1-dimensional 
quantum rotor model can be found in FigfTSl where we made use of finite- iV MPS algorithms with d = 14 and \ = 16. 
The fit to Eq. (|42|) is good and provides a value of c = 0.95 ± 0.09, which is to be compared with the exact known 
value c = 1, corresponding to the £/(l)-invariant critical point. 



We have also performed a numerical evaluation in the thermodynamic limit of the spectrum of the reduced density 
matrix of half an infinite system. First, the convergence of the spectrum for finite x = 30 and off-critical UJ J has 
been considered as the local dimension d increases. Our results are shown in Fig ]13[ where it is possible to observe 
the convergence in the relative error of each one of the x eigenvalues (Aq.) 2 as d increases. This error is found to be 



B. Finite-size corrections to the critical ground state energy 




(42) 



C. Spectrum of half an infinite system 
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of the order of 10 -8 . Second, we have performed a computation of the behaviour of the spectrum of the half-infinite 
reduced density matrix as a function of U/J and away from criticality. In Fig |14l we show the found results for d = 7 
and x = 20, where a remarkable monotonicity in the eigenvalues and their degeneracies when flowing away from 
criticality in the parameter space is appreciated. Indeed, it is possible to prove that such a monotonicity is very much 
related to the preservation of part of the conformal structure of the critical point, when flowing away from criticality 
in the space of parameters, involving in turn the monotonic majorization of the considered probability distribution 

His sum. 



VI. NUMERICAL SOLUTIONS OF PARTIAL DIFFERENTIAL EQUATIONS 

The fact that a set of coefficients can be approximated up to some accuracy by means of the contraction of a given 
tensor network (such as MPS), has applications well beyond the description of quantum many-body systems, as is the 
case of the MPS-assisted image compression 38]. In this section we show that the techniques previously introduced 
can also be used to deal with a purely mathematical problem, namely, that of finding numerical solutions to partial 
differential equations. 

Let $(xi , . . . , xn) denote a function of N real variables defined on some domain D. Let O denote a linear differential 
operator, that is we assume that O is some polynomial in xt,... : xn,9 Xi , ■ ■ ■ ,d XN . The problem that we consider is 
to find a numerical solution, \&(xi, . . . , xn), of the equation 

6* + $ = 0, (43) 

on D. We suggest that for wide classes of partial differential equations, using iV-variable matrix product functions 
might turn a powerful computational method. The intuition underlying our point is most simple. A celebrated 
method to solve partial differential equations is to assume separation of variables, that is, that the solution has the 
form n^L^W (xi). In the language of quantum information, such a solution would correspond to a pure product state. 
On another hand, a fully general solution would have the form 



oo 

. . . ,X N ) = ^2 c ( s lT--' S N)<t>ki( s l)--- ( f> i s^H x N)- 



(44) 
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-0.95 
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-1.15 



FIG. 10: Convergence with N of the energy per lattice site of the ground state for the 1-dimensional quantum rotor model, for 
d — 7, x = 20 and U/J = 1, as computed with MPS algorithms for finite N and in the thermodynamic limit. 



1.4 




FIG. 11: Entanglement entropy of half a system for the 1-dimensional quantum rotor model, for d — 7 and \ — 20, as computed 
with MPS algorithms for finite N and in the thermodynamic limit. Spontaneous symmetry breaking of the Z2 symmetry of the 
ground state is observed for the algorithms with finite N towards the limit U/J^0, which is not observed in the infinite-size 
case. This makes the wave function to decay from a GHZ-like superposition of two states to a completely separable state, which 
results in a null value of the entanglement entropy. 
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FIG. 12: Finite size scaling of the ground state energy at the point of maximum correlations for the 1-dimensional quantum 
rotor model, as computed with MPS algorithms with d — 14 and \ = 16- The line represents the fit to Eq. (|42|l . which provides 
a value of the central charge of c = 0.95 ± 0.09, which is compatible with the exact value c = 1. 



where the tensor c(s±, . . . ,sn) is a priori not expected to have any particular structure. In the language of quan- 
tum information theory, Eq. (|44p would correspond to an arbitrary pure entangled state, and could be, in principle, 
exponentially hard to describe as TV grows. Our point is that, in between, the finitely correlated MPS offer a much 
more refined description of the solution than the mere separation of variables, but still allow for efficient numerical 
computations. 



A. The algorithm 



Let us work in what follows with periodic boundary conditions. In order to find a solution to Eq. (|43p . we will seek 
at minimising the error 



W= I \CW + <f>\ 2 . (45) 
D 



Now observe that if one sees the components of the matrices A(k,Sk) related to the variable k as those of a 
d\ -component vector, a&, then W appears as a positive quadratic form: 



W = &* k ■ R^ ■ a k + r^* • a k + a* • T« + / 



*\, (46) 



where is a dx 2 x d\ 2 positive matrix and r'" is a d\ 2 vector. Our algorithm, inspired from those used to find 
ground states of Hamiltonians described in the previous sections reads as follows: 

• Set the matrices A(k, Sk) to some initial random values. 

• Sweep through each variable, Xk, sequentially, and improve the matrices A(k, Sk) upon solving the linear system 
of equations 

R {k) ■ a k + T (fe) = 0, (47) 
until some desired convergence is attained. 
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FIG. 13: Spectrum of the reduced density matrix of half an infinite system for the 1-dimensional quantum rotor model, as 
computed with MPS algorithms in the thermodynamic limit and \ = 30. In the inset, convergence is seen to appear in the 
relative error of each one of the eigenvalues when increasing the local dimension d. 




a 



FIG. 14: Spectrum of the reduced density matrix of half an infinite system for the 1-dimensional quantum rotor model, as 
computed with MPS algorithms in the thermodynamic limit with d — 7 and \ — 20, when flowing away from criticality in the 
ratio U/J. A remarkable monotonicity in the spectrum is observed. 
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At this point, it is worth informing the reader of three useful tricks that help to make the computations well- 
behaved. First, the variables should always be "visited" in the same order. For example, if the labels of the variables 
were arranged on a ring, the matrices related to each variable should always be improved in a clockwise order. Second, 
we have chosen the random matrices to which the MPS is initialised to be unitary (preconditioning of the solution). 
Third, after improvement of the matrices related to the k-th variable, we recommend to perform the gauge trans- 
formation described in j3]: we compute the singular decomposition A(k, Sfe) Qfc _ l!afc = ^(sla k l ) A^A^^iv** 
where f/Wtf/W = I X) and make the gauge transformation A(k, Sk) ak _ liak — > U(k, Sk) ak - U \ = ^v* a ) \- 



B. An example 

We have considered the Poisson equation: 
A* + $ = 0, (48) 

where A = YliLi ~§^?- The domain we have considered is D = [0, l]i x • • • x [0, l]jy. As individual set of functions, we 
have chosen the plane waves with integer wave index for all the possible sites: 

4>i k J(xk) = <f>s k (x k ) =e i2 ™* x \ s k = -(d-l)/2,...,(d-l)/2 Vfc. (49) 
We choose 3> to be of the MPS form. 

(d-l)/2 

$(x u ...,x N )= tx{F(l,s 1 )...F(N,s N ))4> sl (x 1 )...ct> SN (x N ), (50) 

sl ...s N =-(d-l)/2 

where the matrices F have size \<s> x X*- Note that taking and d large enough, any (smooth) function $ can be 
represented in the form (|50p. 

We have performed three different sets of computations for N — 8, 20 and 40 variables. All our computations 
are such that x$ = d, and such that the matrices F are chosen to be random and unitary. To fairly evaluate the 
performance of our algorithm, we define a relative error 5 as 

+ (51) 



/„ |0*P + /„ |4f + 2^/„ |0*|= x / D |*|» 

This definition is inspired by the Cauchy-Schwarz inequality (0 < 5 < 1). 

Our results are displayed in table [TTJ It is remarkable that the relative error decreases approximately exponentially 
with x- 



C. Including boundary conditions 

We have shown how to solve an equation on a given domain without extra constraints. Therefore, the algorithm 
converges to some solution. Still, one often demands to be able to take into account some boundary conditions. 
Suppose that they are given by a set of m relations of the form 

e i [V(x 1 ,...,x N )} + fi(x u ..., x N ) = (52) 

where the relation i is expected to hold on some domain T>i, i — 1, ...,m. When is again of the form 
polypi, . . . , a; 7v, d Xl , . . . , d XN ), then the boundary condition error 

W BC = y I \e z [y(x 1 ,... 1 x N )] + f t (x 1 ,...,x N )\ 2 (53) 

is of the form (|46p for some matrices R^q, and some vectors rj^,. So after making the transformation — > 

R^ + RftQ, r( fe ) — > r( fc ) + Tg^, we see that, with our method, solving a partial differential equation with boundary 
conditions can be done in exactly the same way as without. 
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N 


d 


X* 


X 


W 


5 


8 


3 


3 


3 


3.25 x 10- 


-3 


8.13 x 10 


-4 


8 


3 


3 


4 


8.4 x 10- 


-4 


2.08 x 10 


-4 


8 


3 


3 


5 


1.36 x 10- 


-4 


3.4 x 10- 


-5 


8 


3 


3 


6 


4.9 x 10- 


-5 


1.23 x 10 


-5 


8 


3 


3 


7 


3.4 x 10- 


-5 


8.51 x 10 


-6 


8 


7 


7 


8 


4.15 x 10- 


-3 


1.04 x 10 


-3 


8 


7 


7 


9 


3.44 x 10- 


-3 


8.61 x 10 


-4 


20 


5 


5 


7 


4.4 x 10- 


-4 


1.1 x 10- 


-4 


20 


7 


7 


3 


4.09 x 10- 


-4 


1.02 x 10 


-4 


20 


7 


7 


4 


3.02 x 10- 


-4 


7.56 x 10 


-5 


20 


7 


7 


5 


1.2 x 10- 


-4 


3 x 10- 


5 


20 


7 


7 


6 


9.73 x 10- 


-7 


2.43 x 10 


-7 


40 


3 


3 


3 


7.3 x 10- 


-5 


1.84 x 10 


-5 


40 


3 


3 


4 


6.17 x 10- 


-5 


1.54 x 10 


-5 


40 


3 


3 


5 


3.96 x 10- 


-5 


9.9 x 10- 


-6 


40 


5 


5 


5 


1.58 x 10- 


-4 


3.94 x 10 


-5 



TABLE II: Performances of the MPS algorithm for solving a Poisson equation. N denotes the number of variables, d denotes 
the number of local levels, x* is the dimension of matrices F, \ is the dimension of the matrices of the MPS ansatz, W is the 
absolute error, and 8 the relative error. 

D. Discussion 

Let us now discuss the principal limitations of the methods. First, when the differential operator and the boundary 
conditions are not linear, the absolute error is no more of the convenient form (|4"fJ)l . More sophisticated methods are 
then necessary to perform the extremisation. 

The second limitation is related to the shape of the domain on which the variables are defined. As a premise to the 
algorithm, one should calculate tensors of the form 

^:::::;£(o.,...,o*)= / 4>*M)i G ^](xi) ■ ■ - ^a^Pn^Kxn), (54) 

where Oi, . . . , On are operators related to the variables 1, . . . , N respectively. 

In the example we have treated, the domain was a mere Cartesian product of sub-domains, which had allowed 
to break G into pieces G^' (Oi, . . . , On) — gil(Oi) • ■ -9°^ (On)- But for an arbitrarily shaped domain, the 
computational resources (time and memory) necessary to compute G could grow exponentially with N. Still the 
method is not limited to the simplest case where the whole domain could be factored as a product of independent 
ones. For example, in the case of problems involving relatively small number of variables, this limitation is not so 
constraining. Second, when D can be decomposed as 

- :i U' :i ,.i - ••• - • ( 55 ) 

i 

where each Djj involves a small number of variables, then G can be decomposed accordingly, and thus computed 
efficiently. These first two comments hold for the boundary conditions as well. 

The other drawbacks of the methods are those generally related to MPS algorithms. First, we have no guarantee 
that the solution provided is a global minimum. However, contrarily to the study of Hamiltonian ground states, we 
have here a figure of merit that tells us how close to the actual solution the algorithm gets. The second drawback 
is related to the very structure of the ansatz. When studying many body quantum systems, it is known that the 
performance of MPS algorithms in their usual context depends on how far one is from criticality, performing well 
when the correlations between the particles are finite. In turn we expect our algorithms to behave similarly, and to 
perform poorly for equations which solution are highly entangled functions. 

Despite its limitations, it is obvious that the method is very versatile and we believe that it can be successfully 
applied to a very wide class of differential problems. It certainly deserves further tests on relevant equations as well 
as theoretical investigation. 
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VII. CONCLUSIONS AND OUTLOOK 



In this paper we have presented in a pedagogical way several numerical optimization techniques based on MPS, 
both for finite systems (with open and periodic boundary conditions) and in the thermodynamic limit. We have 
shown how these algorithms can be extended in order to study quantum many-body systems with continuous degrees 
of freedom (continuous variables), by means of appropriate truncations in the local Hilbert spaces. On the one hand, 
and as a test, we have applied our algorithms to study a finite 1-dimensional chain of harmonic oscillators. On the 
other hand we have performed a detailed numerical study of the properties of the 1-dimensional quantum rotor model. 
In both cases, we have found high precision results when working off criticality. 

The methods presented in our work can be further applied to the study of different continuous-variable quantum 
many-body systems. For example, it could be possible to consider the behaviour of anyonic excitations in small Bose- 
Einstein condensates [43| , by means of a numerical treatment of Laughlin and Pfafnan wavefunctions |4f| 0, [47J . 
Also, it could be possible to perform a study of systems which alternate both discrete and continuous degrees of 
freedom, such as spin-boson models ^48j. Interacting quantum field theories such as the (1 + l)-dimensional \ip 4 
theory can also be studied by means of these techniques [23j], and will be the main topic of a further work [49| . 

Clearly, our work can also be extended to alternative tensor networks, such as the 2-dimensional PEPS [8| or the 
MERAs [13]. The most delicate point when applying MPS algorithms is probably the choice of a basis. In the case 
we have treated, the choices (Fock states or plane waves) came quite naturally. But more often than not, to guess a 
correct single-particle basis in which truncate the Hamiltonian, is not an easy business. 

Finally, we have exhibited an MPS-like method to solve partial differential equations. The results we have got with 
the example treated are quite promising, and we hope that the method will find many applications. Further studies 
of it are currently under investigation. 
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